Introduction au Reservoir Computing avec ReservoirPy¶

Nathan Trouvain
Inria - Mnemosyne




AI4I - 17 janvier 2023

Sommaire¶

  • Concepts et fonctionnalités clés
  • Prédictions de séries temporelles
  • Déployer les capacités génératives
  • Apprentissage en ligne
  • Comprendre les hyperparamètres
  • Démonstration d'application - chute de robot
  • Démonstration d'application - chants de canaris
  • Démonstration d'application - analyse de sentiment

Concepts et fonctionnalités clés ¶

  • Numpy, Scipy, et seulement Numpy et Scipy
  • Exécution efficace (parallélisation des calculs, optimisation des règles d'apprentissage)
  • Plusieurs méthodes d'apprentissage, online, offline
  • Flexibilité des architectures
  • Outils d'aide à l'optimisation des paramètres
  • Documentation: https://reservoirpy.readthedocs.io/en/latest/
  • GitHub: https://github.com/reservoirpy/reservoirpy

Informations générales¶

  • Toutes les données sont des tableaux NumPy
  • Le temps est toujours représenté par la première dimension des tableaux/vecteurs.

Prédiction de séries temporelles ¶

Un exemple: l'attracteur de Lorenz

$$ \begin{split} \dot{x}(t) &= \sigma (y(t) - x(t)) \\ \dot{y}(t) &= \rho x(t) - y(t) - x(t)z(t) \\ \dot{z}(t) &= x(t)y(t) - \beta z(t) \end{split} $$
  • décrit l'évolution des phénomènes de convection dans un fluide
In [2]:
from reservoirpy.datasets import lorenz

timesteps = 3000
X = lorenz(timesteps, x0=[17.67, 12.93, 43.91])
In [4]:
plot_lorenz(X, 1000)

Connaissant la valeur de la séries à l'instant $t$:

  • Comment prédire $t+1$, $t+100$...?
  • Comment prédire $t+1$, $t+2$, $\dots$, $t+n$ ?

Prédire la valeur de la série à un horizon de 10 pas de temps¶

Prédire $P(t + 10)$ connaissant $P(t)$.

In [6]:
from reservoirpy.datasets import to_forecasting

x, y = to_forecasting(X, forecast=10)
X_train1, y_train1 = x[:2000], y[:2000]
X_test1, y_test1 = x[2000:], y[2000:]

plot_train_test(X_train1, y_train1, X_test1, y_test1)

Rappels sur le Reservoir Computing¶

$$ x[t+1]= (1 - \alpha) x[t] + \alpha f(W \cdot x[t] + W_{in} \cdot u[t] + W_{fb} \cdot y[t]) $$$$ y[t]= W_{out}^{\intercal} x[t] $$

Echo State Networks avec ReservoirPy¶

Utilisation de Nodes pour constuire des modèles:

Préparation de l'ESN¶

In [8]:
units = 100               # - nombre de neurones
leak_rate = 0.3           # - leaking rate
spectral_radius = 0.95    # - rayon spectral
input_scaling = 0.5       # - facteur de mise à l'échelle des entrées (input scaling)
connectivity = 0.1        # - densité des connexions du reservoir vers lui même
input_connectivity = 0.2  # et des entrées vers le reservoir
regularization = 1e-4     # - coefficient de régularisation (L2)
transient = 100           # - temps de chauffe du reservoir
seed = 1234               # - reproductibilité
In [10]:
from reservoirpy.nodes import Reservoir, Ridge

reservoir = Reservoir(units, input_scaling=input_scaling, sr=spectral_radius,
                      lr=leak_rate, rc_connectivity=connectivity,
                      input_connectivity=input_connectivity, seed=seed)

readout   = Ridge(ridge=regularization)
In [11]:
esn = reservoir >> readout
In [12]:
reservoir_fb = reservoir << readout

esn_fb = reservoir_fb >> readout

Entraînement de l'ESN¶

L'apprentissage est offline ("hors-ligne") : il n'a lieu qu'une seule fois, sur l'ensemble des données d'entraînement.

In [13]:
esn.fit(X_train1, y_train1, warmup=transient);
In [15]:
plot_readout(readout)

Test de l'ESN¶

In [16]:
y_pred1 = esn.run(X_test1)
In [17]:
plot_results(y_pred1, y_test1)

Coefficient de détermination $R^2$ et erreur quadratique normalisée :

In [18]:
rsquare(y_test1, y_pred1), nrmse(y_test1, y_pred1)
Out[18]:
(0.9961759934656298, 0.012317846485056738)

Déployer les capacité génératives ¶

  • Entraînement de l'ESN sur une tâche de prédiction de 1 pas de temps.
  • Prédiction de l'ESN sur ses prédictions précédentes (mode génératif).
In [19]:
units = 300               # - nombre de neurones
leak_rate = 0.3           # - leaking rate
spectral_radius = 1.25    # - rayon spectral
input_scaling = 0.1       # - facteur de mise à l'échelle des entrées (input scaling)
connectivity = 0.1        # - densité des connexions du reservoir vers lui même
input_connectivity = 0.2  # et des entrées vers le reservoir
regularization = 1e-4     # - coefficient de régularisation (L2)
transient = 100           # - temps de chauffe du reseroivr
seed = 1234               # - reproductibilité

Entraînement à la prévision sur un horizon court¶

In [21]:
esn = reset_esn()

x, y = to_forecasting(X, forecast=1)
X_train3, y_train3 = x[:2000], y[:2000]
X_test3, y_test3 = x[2000:], y[2000:]

esn = esn.fit(X_train3, y_train3, warmup=transient)

Génération¶

  • 100 pas de temps de la série de test utilisés comme "échauffement";
  • 300 pas de temps générés "à partir de rien";
In [22]:
seed_timesteps = 100

warming_inputs = X_test3[:seed_timesteps]

warming_out = esn.run(warming_inputs)  # échauffement
In [23]:
nb_generations = 500

X_gen = np.zeros((nb_generations, 3))
y = warming_out[-1]
for t in range(nb_generations):  # génération
    y = esn(y)
    X_gen[t, :] = y
In [24]:
X_t = X_test3[seed_timesteps: nb_generations+seed_timesteps]
plot_generation(X_gen, X_t, warming_out=warming_out, warming_inputs=warming_inputs)
In [25]:
plot_attractors(X_gen, X_t, warming_inputs, warming_out)

Apprentissage online ¶

Apprentissage se déroulant de manière incrémentale.

Utilisation de l'algorithme Recursive Least Squares

In [27]:
from reservoirpy.nodes import RLS

reservoir = Reservoir(units, input_scaling=input_scaling, sr=spectral_radius,
                      lr=leak_rate, rc_connectivity=connectivity,
                      input_connectivity=input_connectivity, seed=seed)

readout   = RLS()


esn_online = reservoir >> readout

Entraînement pas à pas¶

In [28]:
outputs_pre = np.zeros(X_train1.shape)
for t, (x, y) in enumerate(zip(X_train1, y_train1)): # pour chaque pas de temps de la série :
    prediction = esn_online.train(np.atleast_2d(x), np.atleast_2d(y))
    outputs_pre[t, :] = prediction
In [29]:
plot_results(outputs_pre, y_train1, sample=100)
In [30]:
plot_results(outputs_pre, y_train1, sample=500)

Entraînement sur une séquence complète¶

In [32]:
esn_online.train(X_train1, y_train1)

pred_online = esn_online.run(X_test1)  # Wout est maintenant figée
In [33]:
plot_results(pred_online, y_test1, sample=500)

Determination coefficient $R^2$ and NRMSE:

In [34]:
rsquare(y_test1, pred_online), nrmse(y_test1, pred_online)
Out[34]:
(0.9949772578256109, 0.014117121497420036)

Comprendre les hyperparamètres et leurs effets: plonger dans le reservoir ¶

In [35]:
units = 300               # - nombre de neurones
leak_rate = 0.3           # - leaking rate
spectral_radius = 1.25    # - rayon spectral
input_scaling = 0.1       # - facteur de mise à l'échelle des entrées (input scaling)
connectivity = 0.1        # - densité des connexions du reservoir vers lui même
input_connectivity = 0.2  # et des entrées vers le reservoir
regularization = 1e-4     # - coefficient de régularisation (L2)
transient = 100           # - temps de chauffe du reseroivr
seed = 1234               # - reproductibilité

1. Le rayon spectral¶

Le rayon spectral est la valeur propre maximale de la matrice des poids du réservoir ($W$).

In [36]:
states = []
radii = [0.1, 1.25, 10.0]
for sr in radii:
    reservoir = Reservoir(units, sr=sr, input_scaling=0.001, lr=leak_rate, rc_connectivity=connectivity,
                         input_connectivity=input_connectivity)

    s = reservoir.run(X_test1[:500])
    states.append(s)
In [37]:
units_nb = 20

plt.figure(figsize=(15, 8))
for i, s in enumerate(states):
    plt.subplot(len(radii)*100+10+i+1)
    plt.plot(s[:, :units_nb], alpha=0.6)
    plt.ylabel(f"$sr={radii[i]}$")
plt.xlabel(f"Activations ({units_nb} neurons)")
plt.show()

  • $-$ rayon spectral $\rightarrow$ dynamiques stables

  • $+$ rayon spectral $\rightarrow$ dynamiques chaotiques

2. Le facteur de mise à l'échelle des entrées (input scaling)¶

Il s'agit d'un coefficient appliqué à $W_{in}$, venant changer l'échelle des données en entrée.

In [38]:
states = []
scalings = [0.00001, 0.001, 1.0]
for iss in scalings:
    reservoir = Reservoir(units, sr=spectral_radius, input_scaling=iss, lr=leak_rate,
                          rc_connectivity=connectivity, input_connectivity=input_connectivity)

    s = reservoir.run(X_test1[:500])
    states.append(s)
In [40]:
units_nb = 20

plt.figure(figsize=(15, 8))
for i, s in enumerate(states):
    plt.subplot(len(scalings)*100+10+i+1)
    plt.plot(s[:, :units_nb], alpha=0.6)
    plt.ylabel(f"$iss={scalings[i]}$")
plt.xlabel(f"Activations ({units_nb} neurons)")
plt.show()

Correlation moyenne des activités des neurones du reservoir avec les entrées :

In [41]:
for i, s in enumerate(states):
    corr = correlation(states[i], X_test1[:500])
    print(f"IS : {scalings[i]}, correlation max : {corr}")

IS : 1e-05, correlation max : 5654.146949848663
IS : 0.001, correlation max : 5658.500729694261
IS : 1.0, correlation max : 5804.2333455595935
  • $+$ input scaling $\rightarrow$ activités corrélées aux données
  • $-$ input scaling $\rightarrow$ activités libres

L'input scaling peut aussi être utilisé pour ajuster l'influence de chaque donnée en entrée.

3. Le taux de fuite de la mémoire (leaking rate)¶

$$ x(t+1) = \underbrace{\color{red}{(1 - \alpha)} x(t)}_{\text{état actuel}} + \underbrace{\color{red}\alpha f(u(t+1), x(t))}_{\text{données suivantes}} $$

avec $\alpha \in [0, 1]$ et:

$$ f(u, x) = \tanh(W_{in} \cdotp u + W \cdotp x) $$
In [42]:
states = []
rates = [0.02, 0.2, 0.9]
for lr in rates:
    reservoir = Reservoir(units, sr=spectral_radius, input_scaling=input_scaling, lr=lr,
                          rc_connectivity=connectivity, input_connectivity=input_connectivity)

    s = reservoir.run(X_test1[:500])
    states.append(s)
In [43]:
units_nb = 20

plt.figure(figsize=(15, 8))
for i, s in enumerate(states):
    plt.subplot(len(rates)*100+10+i+1)
    plt.plot(s[:, :units_nb] + 2*i)
    plt.ylabel(f"$lr={rates[i]}$")
plt.xlabel(f"States ({units_nb} neurons)")
plt.show()
  • $+$ leaking rate $\rightarrow$ faible inertie, faible mémorisation des états précédents
  • $-$ leaking rate $\rightarrow$ forte inertie, grande mémorisation des états précédents

Le leaking rate contrôle la "mémoire" de l'ESN. Il peut être vu comme l'inverse de sa constante de temps

Application réelle : chute de robot ¶

In [45]:
features = ['com_x', 'com_y', 'com_z', 'trunk_pitch', 'trunk_roll', 'left_x', 'left_y',
            'right_x', 'right_y', 'left_ankle_pitch', 'left_ankle_roll', 'left_hip_pitch',
            'left_hip_roll', 'left_hip_yaw', 'left_knee', 'right_ankle_pitch',
            'right_ankle_roll', 'right_hip_pitch', 'right_hip_roll',
            'right_hip_yaw', 'right_knee']

prediction = ['fallen']
force = ['force_orientation', 'force_magnitude']
In [50]:
plot_robot(Y, Y_train, F)

Entraînement de l'ESN¶

Utilisation de l'objet ESN, modèle optimisé et distribué pour les Echo State Networks: permet la parallelisation des calculs

In [52]:
from reservoirpy.nodes import ESN

reservoir = Reservoir(300, lr=0.5, sr=0.99, input_bias=False)
readout   = Ridge(ridge=1e-3)
esn = ESN(reservoir=reservoir, readout=readout, workers=-1)  # version distribuée
In [53]:
esn = esn.fit(X_train, y_train)
In [54]:
res = esn.run(X_test)
In [57]:
plot_robot_results(y_test, res)
In [58]:
print("RMSE moyenne :", f"{np.mean(scores):.4f}", "±", f"{np.std(scores):.5f}")
print("RMSE moyenne (avec seuil) :", f"{np.mean(filt_scores):.4f}", "±", f"{np.std(filt_scores):.5f}")

RMSE moyenne : 0.1693 ± 0.10344
RMSE moyenne (avec seuil) : 0.1443 ± 0.15187

Application réelle : le chant des canaris domestiques ¶

Les données peuvent être téléchargées sur Zenodo : https://zenodo.org/record/4736597

Plusieurs motifs temporels répétitifs differents à décoder : les phrases.

  • Un label par type de phrase à classifier dans le temps.
  • Un label SIL annotant le silence, à détecter également pour permettre la segmentation du chant.
In [59]:
im = plt.imread("./static/canary_outputs.png")
plt.figure(figsize=(15, 15)); plt.imshow(im); plt.axis('off'); plt.show()

Entraînement de l'ESN¶

In [65]:
esn = esn.fit(X_train, y_train)
In [66]:
outputs = esn.run(X_test)
In [68]:
scores  # pour chaque chant testé
Out[68]:
[0.9494949494949495,
 0.9085303186022611,
 0.9531759739013625,
 0.9410002304678498,
 0.9638157894736842,
 0.9478260869565217,
 0.892835458409229,
 0.9403111739745403,
 0.9212787899621864,
 0.8765393116514051]
In [69]:
print("Précision moyenne :", f"{np.mean(scores):.4f}", "±", f"{np.std(scores):.5f}")

Précision moyenne : 0.9295 ± 0.02716

Application réelle: analyse de sentiment ¶

Analyse des données de Maas et al. (2011): IMDB dataset.

Composé d'avis attribués à des films par les utilisateurs de IMDB.

Objectif: attribuer un score de 0 ou 1, aux commentaires respectivements positifs et négatifs.

  • Technique proposée: RC sur word embeddings

Utiliser spacy¶

In [70]:
import spacy

nlp = spacy.load("en_core_web_lg")
In [71]:
sad_doc = nlp("I hate everything about this movie.")
happy_doc = nlp("I really liked this movie!")
unrelated_doc = nlp("Dogs usually hate cats.")
In [72]:
print("Similarités: ", sad_doc.similarity(happy_doc), sad_doc.similarity(unrelated_doc))

Similarités:  0.9271912733659942 0.5860950057882411
In [73]:
for token in sad_doc:
    print(token)

I
hate
everything
about
this
movie
.
In [74]:
print(token.vector)

[-0.076454 -4.6896   -4.0431   -3.4333   11.758     3.7212   -0.98133
  2.7902    0.43608  -2.4425    6.9311    2.7647   -1.7035    6.3123
  2.6301    1.1136   -3.3689   -5.4062   -3.6021   -9.8479    4.6016
 -0.50542   0.66145  -3.552    -4.3909   -2.3702   -1.8251   -4.8768
  2.7245    5.3177   -0.50146  -7.2556   -4.6895   -5.4604   -5.6257
  1.6121   -0.8789    1.6936    1.9308    6.7769    3.6739    0.28693
 -0.22849  -4.0779   -4.7086    2.0722   -2.3405   -1.6884   -1.9876
 -2.2207   -2.6514    3.3686    0.12579  -5.8689    4.6715    1.1999
 -4.3008    2.1239    0.029298  3.0072   -2.4435   -3.9103    1.0187
 -2.0263    5.0496    0.66619  -3.978    -4.0066   -1.8538    0.65539
 -2.1525   -9.2576   -2.6783   -1.4023    4.0334   -2.3129   -5.1051
 -0.91791   0.11236   2.626    -5.5713   -2.7705   -2.5144    0.069018
  6.1091   -3.0522    0.025957 -3.9766   -0.42248   0.60821  -4.883
  2.8595    0.046312 -5.0708   -0.87191  -3.0453   -1.5862   -2.7432
  3.6893    0.69812   6.3418    2.264     2.7177    4.6719    2.9297
  1.1045    2.3091   -0.52648  -2.776    -2.0394   -4.1364    0.091809
 -0.62613  -0.14244   1.4294    5.0263    4.395    -0.63146   3.0023
 -3.5179   -0.71451  -7.009    -0.08824   3.5531    3.1061   -6.3624
  6.6491   -4.8624    0.20515   0.17462  -3.203     1.4741    4.2968
 -4.4224    0.31291  -0.34638  -6.8143   -6.1415    3.9493   -5.1503
  0.29465  -2.0811   -2.1067    2.1568   -3.5174    1.6831   -2.5722
 -4.3133    3.0666   -1.7038   -2.353     0.81963   2.6208    5.6669
 -1.8992    0.57382   5.4505    3.6633   -4.7065    0.60155  -0.53312
 -3.1413   -2.4749    8.7393   -1.3324   -4.1015   -2.9628   -0.079662
 -1.5063    3.3633    7.0645    1.2196    4.3788   -6.4489    1.2194
  0.37559  -1.2663    2.947    -2.7332   -3.5665   -1.3088    0.19732
 -0.77669  -6.1613   -4.2444   -0.68401  -0.61397   0.030807 -0.12418
  2.8148   -2.2854   -2.83     -1.9922   -7.3996   -2.8544   -3.7009
 -0.47745   1.5749   -1.9741   -1.5416    0.45149  -0.21663   1.1495
 -5.0475   -0.19798   0.19929  -3.7595    2.5563   -3.8548   -2.7141
  2.9688   -1.3176   -1.2686   -1.3209   -2.7947    3.6765   -6.7201
 -3.5003   -0.1796    1.1076   -1.0369   -0.27029   1.3922    0.1615
 -1.952    -3.9036    1.1277    2.3389    5.1609    2.1905    2.9146
 -8.1347    0.68726   2.8535    1.4085    5.4403   -9.4458    2.8775
 -0.66906  -1.1599   -3.3115   -3.5608    1.5127   -3.1845   -9.1884
  0.6305   -5.6819   -0.12861  -2.3537    1.0302   -0.76796  -0.83891
 -4.1329    1.773    -0.073686 -1.2361    5.0983   -3.4001    0.42573
 -1.9998    1.439     4.7645    6.0452    4.9381    2.3885   -1.0576
 -2.0167   -0.4184   -1.8464    2.4374    2.3647    1.2872   -1.7353
  1.6065   -0.98225   2.762     1.4601   -0.5405    0.26084   8.2027
  0.97069   3.2362   -5.9116    3.7992   -0.31084  -1.3803   -3.0381
  0.19792  -1.6612    0.39299  -3.7749    3.0147    1.8615   -3.9556
  1.134     3.0243   -3.0984    1.304    -0.52699  -1.3622  ]

Taille du jeu de données:

  • train: 1000 avis postifis + 1000 avis négatifs, de longueur variable.
  • test: idem.

Chaque avis a une longueur N, son nombre de tokens. Chaque token est encodé par un vecteur de 300 dimensions.

Entraînement de l'ESN¶

In [79]:
from reservoirpy.nodes import Reservoir, Ridge

reservoir = Reservoir(500, sr=0.9, input_scaling=0.01, lr=0.1, input_bias=False)
readout = Ridge(ridge=1e-8)

esn = reservoir >> readout
In [80]:
esn.fit(X_train, y_train)
Out[80]:
'Model-5': Model('Reservoir-16', 'Ridge-4')

Prédiction¶

On obtient un score par pas de temps. Pour réduire: on prend le maximum de score sur tous les pas de temps, et on applique un seuil.

In [81]:
y_pred = esn.run(X_test)
In [82]:
threshold = 0.85  # choisi grâce à la ROC
In [84]:
plot_esn_response()
In [85]:
y_pred_thr = [1 if y.max() > threshold else 0 for y in y_pred]
y_test_thr = [1 if y.max() > threshold else 0 for y in y_test]
In [86]:
from sklearn.metrics import accuracy_score

score = accuracy_score(y_test_thr, y_pred_thr)
In [87]:
print("Précision: ", score)

Précision:  0.7685

Merci de votre attention.¶

Nathan Trouvain
Inria - Mnemosyne




AI4I - 20 janvier 2021